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1. Introduction 


This project investigates computational modeling of fatigue crack growth in spiral 
bevel gears. Current work is a continuation of the previous efforts made to use the 
Boundary Element Method (BEM) to simulate tooth-bending fatigue failure in spiral 
bevel gears. This report summarizes new results predicting crack trajectory and fatigue 
life for a spiral bevel pinion using the Finite Element Method (FEM). 

Predicting crack trajectories is important in determining the failure mode of a 
gear. Cracks propagating through the rim may result in catastrophic failure, whereas the 
gear may remain intact if one tooth fails and this may allow for early detection of failure. 
Being able to predict crack trajectories is insightful for the designer. However, predicting 
growth of three-dimensional arbitrary cracks is complicated due to the difficulty of 
creating three-dimensional models, the computing power required, and absence of closed- 
form solutions of the problem. 

Previously in this project, tooth-bending fatigue failure in spiral bevel gears was 
investigated using the BEM [1]. These analyses were significant in developing a method 
for predicting three-dimensional, non-proportional, fatigue crack growth incorporating 
moving loads. Prior to the BEM study, there were not many examples of three- 
dimensional crack growth work on gears. Most of the gear studies in the literature are 
two-dimensional analyses. Furthermore, moving loads had not been considered in these 
analyses. 

The BEM predictions were helpful to determine areas in need of further research. 
These were: improving accuracy of simulation by switching from the boundary element 
to the finite element method, decreasing computation time by employing parallel FEM 
analysis, and improving modeling of contact tooth loads. 

Substantial improvement in computation time has been achieved by using parallel 
FEM technologies on PC-clusters. This makes it possible to simulate crack growth in 
large models, like the current gear problem, in a few minutes. It is also possible to obtain 
accurate values of stress intensity factors (SIF) using the FEM and elastic, equivalent 
domain J-integral approaches. A new technique was developed for modeling moving 
tooth loads, which decreases the computation time while introducing a more efficient 
procedure for applying contact loads. 

Another focus of this project was performing three-dimensional contact analysis 
of a spiral bevel gear set incorporating cracks. These analyses were significant in 
determining the influence of change of tooth flexibility due to crack growth on the 
magnitude and location of contact loads. This is an important concern since change in 
contact loads might lead to differences in SIFs and therefore result in alteration of the 
crack trajectory. 

Contact analyses performed in this report showed the expected trend of 
decreasing tooth loads carried by the cracked tooth with increasing crack length. 
Decrease in tooth loads lead to differences between SIFs extracted from finite element 
contact analysis and finite element analysis w'ith Hertz contact loads. This effect became 
more pronounced as the crack grew. 
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2. Three-dimensional FEM Simulations of Fatigue Crack Growth in Spiral Bevel 
Gears 

The main objective of the current work is to demonstrate that we can improve 
accuracy and efficiency of crack growth simulations in spiral bevel pinion through FEM 
computations along with improved modeling of tooth loads. In the current study, as in 
previous BEM simulations, LEFM theories are used for fatigue crack growth predictions 
in gears. Fatigue crack growth rates were determined using a modified Paris model 
accounting for crack closure. Crack trajectory predictions are made using the maximum 
circumferential stress theory. Moving loads were taken into account by discretizing the 
loading into 1 5 steps that created a non-proportional load history at the tooth root. 

Accuracy of computations can be improved using the FEM that allows a better 
SIF calculation method, the equivalent domain J-integral method. Solution time for each 
load step can be decreased substantially by employing parallel FEM analysis. Another 
improvement in the accuracy of the results can be obtained by developing a better 
representation of moving contact loads. This is accomplished by approximating the 
moving contact ellipse by interpolation using the shape functions on the surfaces of the 
loaded elements. In the earlier boundary element computations, contact ellipses were a 
part of the geometry models. Since the ellipses overlapped, four different geometry 
models were needed for each crack configuration. In the current approach, contact loads 
are no longer a part of the geometry model. This approach brings a substantial decrease in 
computation time since only one geometry model for each crack configuration is now 
needed. 

The method used in FEM simulations of spiral bevel gear is composed of the 
following steps: 

1 . Create initial geometry model of the spiral bevel pinion using OSM. 

2. Specify boundary conditions and material properties in FRANC3D. 

3. Initiate crack in FRANC3D. 

4. Create a surface mesh composed of triangular elements in FRANC3D. 

5. Create a 3D FEM mesh of the model composed of tetrahedra by J-Mesh using a 
geometry description file written out from FRANC3D. 

6. Calculate the magnitude and location of the nodal contact loads for each load step 
on the loaded elements using the mesh file produced by J-Mesh. 

7. Run FEM analysis on a parallel PC-cluster for 15 discrete load steps. Using the 
displacement results, calculate SIF values at each discrete crack front point for 
each load step. 

8. Calculate the amount of crack extension at each discrete crack front point as a 
result of 15 steps of non-proportional moving load. 

9. Determine new crack front by piecewise least squares fit of the propagated points 
corresponding to each discrete crack front point in FRANC3D. 

10. Remesh surface locally and repeat steps 5 to 10. 
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2.1 Modeling of Moving Tooth Loads 


Contact in a spiral bevel gear occurs in three-dimensions following a path along 
the tooth surface starting from the fillet of the toe to the top of the heel. In analyzing the 
spiral bevel gear, it is assumed that the contact between the gear and pinion follows Hertz 
theory of elastic contact. Hertzian contact holds as long as the significant dimensions of 
the contact area are small compared to the dimensions of each body and to the relative 
radii of curvature of surfaces [2]. Hertz theory assumes that surfaces are continuous and 
non-conforming, strains are small, each solid can be considered as an elastic half-space 
and surfaces are frictionless. 

Under these assumptions, only normal pressure acts between two bodies due to 
Hertz contact producing normal displacements of the surfaces. Since the contact area is 
assumed to be elliptical the region is defined by, 
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Total load acting on the ellipse and average load are then calculated as, 
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In the above equations p 0 denotes the Hertz stress which is the maximum stress, a and b 
are half of major and minor axes of the ellipse, respectively, Figure 2. 1 . 



Figure 2. 1 : Representation of Hertz stress over an elliptical area. 

In order to determine the contact points on the gear tooth, the method described in 
[3] is used to map the three-dimensional contact path into a two dimensional space. The 
following equations were used for the transformation: 
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Xm,Ym,Zm = coordinates of mean contact point (corresponds to the contact ellipse centers) 
Xn, Yn.Zn — coordinates of an arbitrary contact point 

Xu, Y /' 2 ,Z/ 2 = coordinates of points corresponding to 1 and 2 in Figure 2.2 in three 
dimensions. 


2b 


Figure 2.2: Sketch of transformed contact ellipse with mean contact point M and an 
arbitrary contact point N. 

An algorithm was developed to find the nodes within the contact ellipses and to 
calculate the applied load on these nodes. This algorithm: 

1 . reads the nodal information from the mesh and geometry files created by 3D FEM 
meshing program, 

2. retrieves the geometry information of the n th contact ellipse, 

3. transforms the nodal points defined in 3D coordinates into contact plane, 

4. checks if the mesh nodes are within the contact area bounded by the n th ellipse 
(/i = 1 ... 15), eg. Figure 2.3, 

5. calculates nodal contact loads if the node is within the contact area, and 

6. writes out the information of nodes found by this search and corresponding loads. 




Figure 2.3: Nodal points where contact loads are applied for load step 1 1. 


2.2 Parallel FEM 

Understanding fracture phenomena in materials is crucial for science and 
engineering applications such as crack propagation in spiral bevel gears. However, crack 
growth simulations are complicated and computationally expensive in realistic structures. 
In order to solve realistic 3D problems, it is necessary to develop fast and accurate 
computational simulation tools. 
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The BEM is an alternative to the FEM for solving large problems. The BEM is 
advantageous since it only requires a surface mesh and has fewer equations to be solved. 
However, it entails a fully dense, nonsymmetric system of equations, and it is difficult to 
extract accurate near-crack-front stresses and displacements. In this respect, it is more 
desirable to use FEM. This brings the challenge of implementing fast and robust parallel 
sparse solvers. 

A parallel PC-cluster in Cornell Theory Center was used for the development of a 
parallel FEM solver for crack growth simulations. The hardware consists of 32 nodes of a 
64 node cluster of Dell PowerEdge 1650 servers with a Giganet interconnect running 
Microsoft Windows 2000 Advanced Server. Each node contains 2 Intel Pentium 111 
processors clocked at 1 GHz and 2 GB of RAM. Only one processor per node is used for 
the simulations. 

Parallel FEM simulation of crack propagation makes it possible to simulate large 
problems that were impossible or impractical to solve previously. The objective in 
developing such a system is to be able to solve a real-life problem with multiple crack 
propagation steps with 1 ,000,000 degrees of freedom in an hour. Such a problem with 
10,000 degrees of freedom takes about 100 hours using a state-of-the-art single processor 
workstation [4]. 

Crack propagation in a spiral bevel gear stands as a large, realistic engineering 
problem that incorporates the computational difficulties mentioned above. In that respect, 
the development of parallel FEM simulations is a significant step in achieving the goal of 
developing a practical and accurate numerical design tool for gears. 

2.3 SIF Calculations 

Complex problems, such as the current gear problem, require computation of 
mixed mode SIFs due to the combined loading that they receive. In the literature, several 
approaches for computing SIFs can be found. Among these methods, the displacement 
correlation method [5-8] and equivalent domain J-integral methods [9-14] are widely 
used in numerical computations. 

In the previous BEM simulations of the spiral bevel gear, the displacement 
correlation method was used for calculating SIFs. This method only requires 
displacement information on the boundary near the crack front and is computationally 
efficient. Since BEM solves for displacement information only on the boundaries, this 
method is suitable for the BEM analysis. The displacement correlation method makes use 
of numerical solutions for displacements in combination with the analytic solutions of 
SIFs. This method is only as accurate as the computed displacements at the correlation 
points. 

The equivalent J-integral method is adaptable to complex problems and presents 
much more flexibility in the type of applications for which it can be used. We used the 
equivalent domain J-integral approach instead of displacement correlation in the FEM 
simulations of the spiral bevel gear. 

J-integral is the energy release expression associated with crack growth defined as 
an integral along a contour that encloses a crack tip. J-integral is defined as [15], 


NASA/CR 2003-212529 


5 



where W is the strain energy density, rij is the unit outward normal to the contour, Gy is 
the traction exerted on the body enclosed by T and crack surface. This integral is 
evaluated along an arbitrary contour enclosing the crack tip and is path independent. 



r r = p + T| + r, + r 


Figure 2.4: Paths of J-integral calculations at a crack tip. 

Numerical evaluation of the J-integral involves the use of equivalent domain 
integral approach. In this approach, the integration along the contour is replaced by an 
integral over a finite size domain. The expression for J-integral then becomes, 


J k = j[ w/8 <,- a , — 

r, c)x k 


k= 1,2. 


( 6 ) 


A continuous function q(x i,X 2 ) is introduced in the above equation in order to replace the 
integration contour T by Tr. Function q has a unit value at the crack tip and it is zero 
along F, Ti and P. The variation of this function within the domain is arbitrary. Using 
divergence theorem equation 6 takes the form. 
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Note that the second term of this integral vanishes for elastic problems. SIFs are 
calculated within the vicinity of the 3D crack front by the pointwise values of J-integral. 
3D cracks can be represented by 2D plane strain fields at the crack front points. These 
fields are satisfied only within close proximity of the crack tip, however it is well known 
that finite element calculations of the field quantities are not very accurate close to the 
crack tip. Equivalent domain integral approach circumvents this difficulty by replacing 
the line integral by a finite domain integral. The relation between the J-integral and SIFs 
is given for an isotropic, linear elastic material as follows, 
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(8) 


J,=^(Kr + K,, 2 ) and J 2 = 

oC j oC/ 

In the above equations J 2 is not path independent. For a pure mode I or mode II problem, 
SIFs can be calculated directly from these equations however under mixed mode loading 
Kj and Ku cannot be calculated separately. Bui [16] proposed a mixed mode separation 
method by decomposing the elastic displacement field into symmetrical and anti- 
symmetrical parts. This results in separation of J into two components that are given as 
follows [16], 


J - J x + J 2 , 






I . W 
U. = U + It; 


(9) 


Examples of numerical investigations on mixed mode problems using this approach 
yielded highly accurate results [17,18]. 

Numerical techniques generally provide good results in an average sense over the 
domain. However, governing equations and/or boundary conditions are only 
approximately satisfied at any one point. This indicates that displacement correlation 
SIFs are not very accurate since they depend on the displacement values at a point. On 
the other hand, J-integral SIFs give better results due to considering much more of the 
results domain over which approximation errors tend to cancel. 

Two models with known analytic solution are investigated to show the higher 
accuracy of the FEM calculations. The first model is a single, internal, flat, penny- 
shaped crack in an infinite body. The radius of the penny-shaped crack is 0. 1 . This crack 
is centered in a block with dimensions 10x5x5. 


t T- 



Figure 2.5: Flat, internal penny-shaped crack subjected to remote tension. 

For this problem Mode II and Mode III SIFs are zero. The analytical solution for 
the Mode I SIF is given by [19], 



(10) 
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Table 2. 1 shows the comparison of solutions from the BEM, and FEM and the analytical 
solution. BEM results shown in Table 2.1 are calculated by the displacement correlation 
method using linear quadrilateral elements with 12 subdivisions along the radius and 
along the quarter of the crack front. FEM calculations are carried out using the equivalent 
domain J-integral method with 36,583 tetrahedra. Note that an elastic modulus of 10,000 
and a Poisson’s ratio of 0.0 are assigned to the model, and a uniform, far- field tension of 
unity is applied. 

Table 2.1: Average Mode 1 SIF values along the crack front calculated by the BEM and 
the FEM compared with the analytical solution. 


Solution method 

Penny-shaped 

% Error 

Analytical 

0.357 

- 

BEM 

0.348 

2.6 

FEM 

0.352 

1.4 


The second example is a single, internal, flat, elliptical crack in an infinite body 
with an elastic modulus of 10,000 and a Poisson’s ratio of 0.0. The orientation of the 
crack is the same as in Figure 2.5. An elliptical crack with 0.1 minor and 0.2 major axis 
dimensions is embedded in a 10x5x5 block. A uniform far-field tension of unity is 
applied. The analytic solution of Kj involves elliptic integrals whereas AT// and Km are 
zero. Detailed discussion of the analytical solution can be found in [19]. Comparison of 
solutions from the BEM, and FEM and the analytical solution are given in Table 2.2. 
BEM calculations are done using the displacement correlation method with linear 
quadrilateral elements with 5 subdivisions along the axes and 12 subdivisions along the 
quarter of the crack front. The FEM results are computed using J-integral method with 
13,464 tetrahedra. 

Table 2.2: Maximum and minimum Mode I SIF values along the crack front calculated 
by the BEM and the FEM compared with the analytical solution. 


Solution method 

Max Ki 

% Error 

Min K, 

% Error 

Analytical 

0.463 

- 

0.327 

- 

BEM 

0.435 

5.9 

0.304 

7.0 

FEM 

0.442 

4.5 

0.333 

1.8 


2.4 Determining New Crack Front 

Crack extension due to one non-proportional load cycle is computed using the 
approach developed in [1]. As previously mentioned, a modified Paris model was used 
that accounts for crack closure effects. A crack is assumed to advance when its SIF is 
large enough to overcome closure and is larger than the SIF of the previous load step. 
The new crack front is determined as follows: 
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1 . 


For every point along the crack front: 

a) Calculate amount of crack growth, da corresponding to each load step i= 
L...15. 


i-l 
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2 . 

3. 


b) Calculate the angle of crack growth corresponding to each load step using 
maximum circumferential stress theory 


0' =2 tan 


I*4±i (A—) 2 + 8 

4 k,: 4 \ k 


c) Calculate the final coordinates of the crack front and trajectory angles by 
approximating the contributions from each load step by a straight line. 

d) Determine number of cycles, N, necessary to have a reasonable amount of 
crack advance compared to pinion’s geometry. 

Update geometry using the new crack front points. 

Repeat 1 and 2 for each crack growth step. 


2.5 Results 


Three-dimensional finite element analyses were performed for an OH-58 spiral 
bevel pinion for 54 crack growth steps, Figure 2.6. In the previous BEM-based analysis, 
only 13 growth steps were performed. One crack growth step consists of 15 static finite 
element analyses in order to simulate moving loads on the pinion tooth. After the 54 th 
crack growth step the analysis was stopped because it was concluded that propagating the 
crack further would not provide additional insights into the simulations and results. All 
predictions were based on the methodology outlined in section 2. 
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Crack 


Figure 2.6: A spiral bevel pinion model with surface FE mesh (a) whole pinion (b) 
up view of pinion teeth. Note the location of the initial crack in the middle tooth. 


Figures 2.7 through Figure 2.9 show Mode 1, II and III stress intensity factors for 
the initial crack configuration for the first eleven load steps. Load steps 1-4 correspond to 
double tooth contact, load steps 5-11 correspond to single tooth contact and load steps 
12-15 again correspond to double tooth contact. Since a crack is not assumed to advance 
when its SIF is smaller than the S1F of the previous load step, only the first 1 1 steps 
practically contribute to crack growth calculations. In these figures crack front position 
corresponds to the points near the crack front at which SIF values are evaluated. In the 
initial step of crack growth simulations the crack front was discretized such that there 
were 222 points at which SIF values were calculated. This is compared to only 61 points 
previously used with the BEM model [1]. 

Figure 2.10 shows the crack trajectory predictions by the FEM on the tooth 
surface and cross section of the tooth for several crack growth steps including the initial 
and final configurations of the crack. Figure 2.1 1 is a close-up view of the last step of 
crack growth from the toe end of the tooth. 



Load step 1 Load step 3 Load step 5 Load step 7 — Load step 9 Load step 1 1 

Figure 2.7: Initial Mode I SIF distribution for load steps one through eleven. Note that 
orientation is from heel to toe. 


NASA/CR 2003-212529 


10 



Figure 2.9: Initial Mode III SIF distribution for load steps one through eleven. Note that 
orientation is from heel to toe. 
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Initial Step: Step 0: 



Step30: 



Midpoint of crack where cross 
section is taken. 






Step 45: 



Figure 2.10: Crack trajectory predictions by FEM on tooth surface and tooth cross section. 


Final Step: Step 54: 
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Figure 2.11: Close-up view of the crack after the last step of crack growth. 

2.6 Comparison of FEM and BEM Results 

This section compares FEM and BEM results in terms of accuracy and efficiency. 
BEM results presented in this section are from reference [1]. One of the most significant 
improvements in the FEM approach over the BEM simulations is the difference in the 
solution time of the problem. This comparison is summarized in Table 2.3. 

Table 2.3: Comparison of performance between FEM and BEM simulations. 


Results 

FEM 

BEM 

Computer power 

32 nodes" 
PC-cluster, 1 GHz 

Single DEC-Alpha 
Workstation, 175 MHz 

Number of unknowns 

-1-3 million 

6700 

(coefficient matrix) 

(sparse, symmetric) 

(dense, unsymmetric) 

Solution time per load cycle 

5-15 hours 

-60 hours 

(1 load cycle=15 load steps) 

(20min- 1 h/load step) 

(4 models) 

Typical time for one crack step 

-8- 1 8 hours* 

-65 hours* 


"Only one processor per node is used. 

Solution time plus 3 hours post-processing and remeshing 
'Solution time plus 5 hours post-processing and remeshing 


The performance comparison in Table 2.3 clearly shows that there is a very 
significant improvement in computation time with the development of a parallel FEM 
solver. The solution time for BEM simulations was about 65 hours for one crack growth 
step whereas FEM simulations took much less time for each crack growth step even 
though the models have many more degrees of freedom. The number of unknowns for the 
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FEM was about a million at the initial step and increased to more than 3 million as the 
crack front advanced. 3D volume meshing of one model ranged from 20 minutes for the 
smallest model to about an hour for the largest models. 

Improvement in the computation time for each crack growth step allowed us to 
increase solution resolution by taking smaller crack increment steps and to simulate more 
steps of propagation. Figures 2.12 to 2.14 show the comparison of crack area, depth and 
crack front length by FEM and BEM. These plots indicate that the FEM-predicted fatigue 
growth rate is marginally lower when we compare the crack area and the crack front 
length and more cycles are required to achieve the same amount of crack growth when 
compared with previous BEM-based results. Figures 2.15 to 2.17 show the SIF 
comparisons between FEM and BEM for load steps 3, 7, and 1 1. 

Table 2.4: Comparison of crack growth amount between FEM and BEM simulations. 


Results 

FEM 

BEM 

Number of steps 

54 

13 

Total number of cycles 

383,000 

311,000 

Final crack front length 

1.55 in 

1 .42 in 

Final crack area 

0.319 in 2 

0.186 in 2 

Final crack depth 

0.268 in 

0.188 in 



Number of cycles 

Figure 2.12: Comparison between crack area versus number of cycles predicted by BEM 
and FEM analyses. 
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0 -I 1 1 1 1 1 1 1 1 

0 50000 100000 150000 200000 250000 300000 350000 400000 

Number of cycles 

Figure 2.13: Comparison between crack depth versus number of cycles predicted by 
BEM and FEM analyses. 



Figure 2.14: Comparison of crack front length between BEM and FEM analyses. 
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18000 


Load Step 3 


13000 A 


r 8000 -l 



-7000 -L 


Crack front position (in) 


—o — Mode I.FEM 

© - Mode I, BEM 

-o- ModeII,FEM 

- o--ModeII,BEM 

— o— Mode III, FEM - 

-© -ModeHI.BEM 


Figure 2. 15: Comparison of initial Mode I, II and III SIFs computed by the FEM and the 
BEM for load step 3. 



— o — Mode I, FEM 

© -Mode I.BEM 

—a— Mode II.FEM 

a - Mode II,BEM 

-°- Mode III,FEM 

© Mode III, BEM 


Figure 2. 16: Comparison of initial Mode I, II and III SIFs computed by the FEM and the 
BEM for load step 7. 
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18000 


Load step 1 1 



-O- Mode I,FEM 

-o- Mode I, BEM 

-o-Mode II, FEM 

"D - Mode II, BEM 

— o— Mode III.FEM 

-o--Mode III,BEM 


Figure 2.17: Comparison of initial Mode I, 11 and 111 SIFs computed by the FEM and the 
BEM for load step 1 1. 



Mode I.FEM o Mode I, BEM — Mode II, FEM 
••a - Mode II, BEM Mode III,FEM ■ o-Mode III,BEM 


Figure 2.18: Comparison of typical S1F history computed by the FEM and the BEM for a 
load cycle at the midpoint of the initial crack front. 
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Figure 2.19: Comparison of typical non-proportional load cycle at the midpoint of the 
initial crack front computed by the FEM and the BEM. 

Figures 2. 1 8 and 2.19 present the comparison of S1F variation over one load cycle 
and non-proportional load history. In Figure 2.18, it is seen that K / predictions are lower 
in the FEM analysis than the BEM analysis over most of the steps. As seen in Figure 
2.19, the ratio of Ku to K/ over 15 load steps shows a similar trend as observed in the 
BEM analysis with a slight difference in values. 



Figure 2.20: Comparison of crack trajectory predicted by the BEM, the FEM and an 
experiment on the tooth surface and tooth cross section at the midpoint of crack. 
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Crack trajectory comparisons in Figure 2.20 show that, on the tooth surface, crack 
trajectory predicted by FEM is similar to the BEM results. The deviation of crack path at 
the toe end in the BEM predictions was also observed in FEM predictions. On the heel 
end, FEM results exhibit a less steep propagation angle than the experimental trajectory. 
On the tooth cross-section it is observed that the predicted crack trajectory by FEM is less 
steep than BEM and experimental results. FEM results also exhibit a shallower ridge 
formation than the experiment. Although the FEM-predicted crack path still does not 
exactly represent the crack trajectory observed in experiments, it captures the global 
behavior of the crack propagation trajectory. 


3. Three-dimensional Finite Element Contact Analysis 

The meshing of a gear and pinion produces loads on the tooth surface that change 
in magnitude and position in time. In section 2.1 analyses assuming Hertzian contact 
theory were described. This assumption, however, is a simplification for a number of 
reasons: 1) the teeth and rims that transmit the load are flexible, 2) assumed elliptical 
contact area on the tooth surface may not be accurate since the curvatures of the teeth 
surfaces in contact are not constant over the contact zone, 3) the center of contact, where 
the maximum contact stress occurs, may differ from the theoretical contact point, and 4) 
most importantly for the current investigation, the tooth flexibility changes as the crack 
propagates in a tooth. The change in flexibility due to cracks may shift the contact area, 
and change the magnitude of the contact loading. 

In this section a series of three-dimensional finite element analyses that model the 
contact conditions between the gear and pinion explicitly are described. The purpose of 
this study is to assess what influence the details of the load transfer through the contact 
conditions have on the computed stress intensity factors, and through them the predicted 
crack trajectory. 

Three-dimensional finite element contact analyses in spiral bevel gears have been 
undertaken by several researchers to study the stresses induced at the root and on the 
surface of the tooth. Earlier work by Bibel et al. [20] utilized gap elements to simulate the 
contact boundary conditions and to evaluate the stresses on the tooth surface. Further 
study on this subject eliminated the need of using gap elements through the use of contact 
capabilities of commercial finite element programs. Preliminary results showing the 
moving contact patches over the tooth surface are reported in Bibel and Handschuh [21]. 
These analyses also incorporate automatic meshing of the gear and pinion via the 
addition of user subroutines to a commercial finite element program. The same authors 
performed additional analyses in order to compare experimentally measured tooth 
bending stresses of spiral bevel gears to FEM results [22]. Other studies on the contact 
analysis of spiral bevel gears, including the development of a methodology combining a 
surface integral and a finite element solution, are described by Vijakaar [23]. Numerical 
examples carried out using this approach show the change in the contact area as a result 
of the misalignment of the gears. 
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3. 1 Approach 


The goal of the present study is to perform contact analyses incorporating cracks 
to determine the effect of changing tooth flexibility on the magnitude and location of the 
contact load and therefore on the stress intensity factors. A commercial finite element 
program, ABAQUS [24], was used to perform the contact analysis in conjunction with the 
software developed by the CFG to calculate fracture-related parameters. A general 
outline of the approach is as follows: 

1 . Create initial geometry model of spiral bevel pinion and gear with OSM using 
the geometry information provided by NASA Glenn Research Center. 

2. Specify boundary conditions, material properties, and contact surfaces in 
FRANC3D. 

3. Create a surface mesh composed of triangular elements in FRANC3D. 

4. Create a 3D FE mesh of the model composed of tetrahedra using the J-mesh 
program. 

5. Transform the FRANC3D model into an ABAQUS input file. 

6. Perform contact analysis in ABAQUS and determine the location and magnitude 
of the contact loads. 

7. Calculate SIFs using the results of FEM contact analysis and also using the results 
from analysis with Hertzian contact loads. Compare these values and investigate 
their effects on the crack trajectory predictions. 

Geometries of a one-tooth sector of the spiral bevel pinion and gear were created 
separately in OSM using the information provided by NASA. This information regarding 
the geometric positions of points defining the gear set was determined analytically using 
the method developed by Litvin et al. [3]. The one-tooth gear and pinion models were 
created in a common coordinate system with the same axis of rotation (z-axis). 
Additional teeth of the gear and pinion were created by copying and rotating the original 
tooth by an angle equal to 360° divided by number of gear teeth. The resulting models 
were further rotated to orient them in the proper meshing position. At the meshing 
orientation there was no interference between the gear and pinion teeth. 

The contact analysis was performed using the contact pair approach in ABAQUS. 
Contact pairs, named as master and slave surfaces, define surfaces which can potentially 
come into contact. If one of the contact surfaces is stiffer or has a coarser mesh, that 
should be selected as the master surface. In this formulation contact direction is always 
normal to the master surface. In ABAQUS, contact can occur in the form of small-sliding 
or finite-sliding. Analyses performed in this work used the finite-sliding formulation 
where surfaces are allowed to undergo finite separation and sliding. This formulation 
defines the master surface as composed of element faces. Continuity of the surface 
normals was attained by smoothing the transitions between the elements. Contact surface 
tractions were calculated assuming that the surfaces are perfectly hard, which means 
there will be no penetration of surfaces into one another. Contact tractions were defined 
by the local basis system formed by the normal to the master surface. 

The analyses used both FRANC3D and ABAQUS. FRANC3D was used to generate 
the complicated geometry of the gear, to specify the analysis attributes, and to define the 
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crack. Also, it was used as a post-processor to calculate the SIFs. ABAQUS was utilized 
for finite element contact analysis of the model. In order to perform contact analysis on 
models with cracks, features that complement FRANC3D and ABAQUS were needed. In 
this respect FRANC3D was extended to include a capability of specifying the master and 
slave contact surfaces and to include a translator to convert the information regarding the 
geometry, mesh, boundary conditions, material properties and contact surfaces to an 
ABAQUS input file. Also, equivalent domain J-integral implementation was added to 
FRANC3D to obtain more accurate SIF values. 

Contact analyses involve nonlinear geometric analysis. The gear/pinion model is 
meshed with linear tetrahedra. Volume meshing is done using J-mesh as before. The 
mesh is finer on the contacting tooth surfaces of both gear and pinion compared to the 
rest of the model. In performing the analyses, it was assumed that there was no friction 
between the contacting surfaces. 

3.2 Analysis 

Finite element contact analysis of a spiral bevel gear and pinion was performed on 
a four-tooth model composed of two pinion teeth and two gear teeth. This model is 
shown in Figure 3. 1 and is adapted from [22]. Only four teeth of the gear and pinion were 
explicitly modeled since this is the minimum number of teeth that is needed to simulate 
double tooth contact. We chose to model only a small portion of the gear and pinion to 
decrease the computation time and effort. This model was adequate to simulate the 
contact conditions and also demonstrate how the stress intensity factors change with 
increasing crack length. 

The model was designed such that load applied at a known distance from the axis 
of rotation transmits the required amount of torque. The gear model was extended to 
include a solid, "wedge shaped" part. This portion was a modeling artifact added so that 
the gear could be allowed to rotate freely about its axis of rotation and come into contact 
with the pinion as a result of the applied torque. As seen in Figure 3.1a, the pinion teeth 
were fixed with zero displacements in x,y,z directions on the rim surfaces. Displacement 
components of the gear were fixed in x,y,z directions along the line of axis of rotation. 
Gear rim surfaces were only constrained in the z direction; x and y degrees of freedom of 
the gear were not fixed in order to allow for free rotation in that plane. Dimensions of the 
gear sector were chosen such that the bending effects caused by the application of load is 
minimized. Early studies showed that when smaller portions of the gear were modeled 
contact loads changed with the change in the location of the applied load due to the 
excessive deformation occurring in the gear hub. 

Full design torque for the pinion was 3099 in-lb. In order to obtain the torque that 
is transferred by the gear this value was multiplied by the gear ratio, giving 1 1580 in-lb. 
A point load was applied at a distance 3.58 in. from the gear axis of rotation to produce 
this amount of torque. The location of the applied load is shown in Figure 3.1a. 
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(b) 



(c) 


Figure 3.1: Contact model of spiral bevel gear set. (a) whole model (b) close up view of 
gear and pinion teeth in mesh (c) FE mesh of the model. Note that two teeth of pinion 
and gear are explicitly modeled. Note also the refined mesh on the pinion and gear teeth 
surfaces where contact is expected to occur. 
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The original configuration of the gear and pinion corresponds to a double tooth 
contact load transfer. In order to obtain a single tooth contact configuration the gear set 
was rotated equivalent to a 12° pinion rotation. During the analyses, first a very low 
displacement was applied to the gear in order to ensure contact between gear and pinion 
and to fix the model in space. Then the load was increased to its full design value. 

In order to obtain reasonable contact loads, an adequately refined mesh was very 
important. Since the gear has a doubly curved geometry, a fine discretization was needed 
to model the curved edges. In our studies, we observed that coarse meshes lead to 
inaccurate contact stresses and stress concentrations. To this end, we discretized the 
surfaces where contact is expected to occur fine enough to obtain reasonable results 
without excessive computational cost. 

3.3 Results 

Contact analyses were carried out on four models which have: no crack, an initial 
crack, and two crack sizes, crack #1 and crack #2. Crack #1 is within the fillet region of 
the pinion tooth. Crack #2 is a larger crack extending to the surface of the tooth. These 
four geometries were analyzed for double tooth contact and single tooth contact 
configurations that correspond to 0° and 1 2° pinion rotation, respectively. 

The number of nodes and elements for the analyzed models ranged between 
12,000-16,000 and 52,000-70,000, respectively. Solution time of a finite element contact 
analysis is about 30 minutes on a 1 GHz Intel Pentium III machine. 




(C) 

Figure 3.2: Pinion teeth with a) initial crack b) crack #1 c) crack #2. Note that in all the 
models, crack is in tooth 1 . 
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3.3.1 Comparison of contact load and area calculated by Hertz theory and three- 
dimensional FEM contact analysis 


Results of contact analyses and how they compare to the Hertz theory are 
summarized in Table 3.1. Total contact load calculated by finite element analysis 
compared well with the total load predicted by Hertz theory. However, distribution of this 
total contact load between two adjacent teeth shows some discrepancies. There is nearly a 
10% difference in the values between contact analysis and Hertz theory predictions on 
each tooth. 

We were unable to attain single tooth contact conditions using this model. At all 
rotations, there is always at least a very little portion carried by the adjacent tooth. This 
might be due to the inaccuracies introduced by the geometric model of the gears. The 
closest case to single tooth loading was obtained when a 12° pinion rotation was 
introduced to the model. This case was assumed to correspond to load step 1 1 during 
which the highest contact load occurs. 

Contact loads calculated by finite element analyses decrease on tooth 1 , where the 
crack is located, and increase on tooth 2 as the crack grows. These results indicate that as 
the crack becomes larger, the flexibility of the pinion tooth is enhanced, therefore the 
contact loads taken up by the two adjacent teeth are redistributed. A larger portion of the 
total contact load has to be carried by the uncracked neighboring tooth. This trend is in 
agreement with the expected behavior. 


Table 3.1: Comparison of contact loads and areas between FEM contact simulations and 
Hertzian stress calculations for two contact configurations with different crack lengths. 
Note that in the first configuration the Hertz contact load given below corresponds to load 
step 3 for tooth 1 and load step 13 for tooth 2, and in the second configuration it 
corresponds to load step 1 1 . 


Model Contact Load (lb) Contact Area (in 2 ) 

Tooth #1 Tooth #2 Tooth #1 Tooth #2 


Configuration # 1 {Rotation 0") 


No crack 

1800 

1399 

0.0129 

0.0126 

Initial crack 

1794 

1404 

0.0129 

0.0126 

Crack #1 

1792 

1406 

0.0123 

0.0126 

Crack #2 

1743 

1457 

0.0177 

0.0126 

Hertz Solution 

1584 

1602 

0.0114 

0.0146 

Con figuration # 2 {Rotation 12" ) 




No crack 

2899 

318 

0.0193 

0.0018 

Initial crack 

2885 

322 

0.0193 

0.0018 

Crack #1 

2881 

326 

0.0193 

0.0018 

Crack #2 

2860 

348 

0.0252 

0.0018 

Hertz Solution 

3196 

- 

0.0216 

- 
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(a) 




(b) 
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(c) 

Figure 3.3: Contact pressures in on pinion tooth surface for a double tooth contact 
configuration with (a) initial crack (b) crack #1 (c) crack #2. Note that the values on the 
contour bars are in psi. 

For earlier steps of the crack growth, when the crack is still in the fillet of the 
pinion tooth, as in crack #1, there is no significant change in the contact load and area on 
the tooth as seen in Table 3.1. However, when the crack extends into the tooth surface, as 
in crack #2, its effect on the flexibility of the tooth is more pronounced. This results in a 
more significant decrease in the contact load on tooth 1. 

There is no significant difference in the centroid of contact loads calculated by 
finite element contact analysis for different crack lengths. However, locations of contact 
loads determined by these analyses were shifted towards the toe compared to the Hertz 
theory predictions. 
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(a) 



(b) 



(c) 


Figure 3.4: Contact pressures on pinion tooth surface for an almost single tooth contact 
configuration with (a) initial crack (b) crack #1 (c) crack #2. Note that the values on the 
contour bars are in psi. 
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Contact areas shown in Figures 3.3 and 3.4 are almost elliptical, showing the 
same trend as the Hertz contact predictions. Also, the distribution of the contact loads is 
close to an ellipsoid shape predicted by Hertzian contact. Maximum contact loads 
calculated by contact analysis were higher than the predictions by the Hertz theory. 

Results of contact analyses are very sensitive to the mesh size of the contacting 
surfaces. If a coarse mesh is used, computed contact pressure is not continuous on the 
surface since only certain nodes are in contact. In order to overcome this issue, the mesh 
on tooth surfaces was refined several times until a continuous contact pressure 
distribution was obtained. The slightly wavy pattern observed in the contact areas in 
Figures 3.3 and 3.4 are an artifact of the mesh size, shape and orientation. Meshes on the 
gear and pinion tooth surfaces are built such that there is compatibility of the mesh nodes 
between these two contact surfaces. This is important since the contact is calculated 
based on the interaction of mesh nodes between contacting surfaces. If the nodes of one 
of the surfaces fall between the two nodes defining an edge of the other surface, then the 
contact load is taken up by either of the nodes that might lead to different contact 
patterns. The discontinuous contact load pattern observed on crack #2 (Figure 3.3c and 
3.4c) for both configurations might be due to the surface mesh of the tooth. At this crack 
step, the crack was extended to the surface of the pinion tooth. This distorted the 
regularity of the mesh on the pinion surface. Since the compatibility of the meshes on 
both surfaces is very important for the contact analysis, the loss of this compatibility 
might lead to the formation of the contact areas as seen in the above figures. However, 
this issue does not bear a high significance for the calculation of the SIFs. For the 
accuracy of the SIF values, the dominating variable is the total contact load transferred 
rather than the exact distribution. 

Figures 3.5 and 3.6 show the comparison of location of contact loads by Hertzian 
analysis and FE contact analysis for both single-tooth and double-tooth contact at initial 
crack. Gray squares in the figures define the extreme points of the Hertzian contact 
ellipses. Black points are the mesh nodes where contact occurs in the FE contact analysis. 
The location of contact loads is almost the same at crack #1 and crack #2. 



Figure 3.5: Comparison of location of contact loads by Hertzian analysis and FE contact 
analysis for double-tooth contact at initial crack. Hertzian loads on tooth 1 and tooth 2 
correspond to load step 3 and load step 13, respectively. 
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Figure 3.6: Comparison of location of contact loads by Hertzian analysis and FE contact 
analysis for double-tooth contact at initial crack. Hertzian load on tooth 1 corresponds to 
load step 1 1 . 

3.3.2 Comparison of SIFs extracted from three-dimensional FEM contact analysis and 
finite element analysis with Hertz contact loads 

Change in contact loads due to cracking is expected to have an influence on the 
stress intensity factors and therefore on the crack trajectory. In order to investigate this 
issue, SIFs were calculated using the results of three-dimensional finite element contact 
analysis and finite element analysis using Hertzian contact loads. The latter set of 
analyses were carried out using a two-tooth pinion model where the Hertzian contact 
loads, corresponding load step 3 and load step 1 1 of the fifteen discretized moving load 
steps described in section 2.1, were applied to the model. 

SIFs were first calculated by the displacement correlation technique for both 
cases; however, contact analysis results did not compare well with the Hertzian contact 
load results. A 30 % difference in SIFs was observed between the two set of analysis. We 
used linear elements for contact calculations, which is required for this type of analysis, 
and extracted SIFs from these results by the displacement correlation technique. This 
approach is not accurate enough with the linear elements. 

The equivalent domain J-integral method is known to give more accurate results, 
and this method was implemented in FRANC3D. However, SIFs obtained by linear 
elements were still somewhat less accurate than the quadratic elements used for the finite 
element analysis with Hertz contact loads. Improved results with the linear elements were 
obtained as the size of the radius determining the domain of the J-integral calculations 
increased. 

Comparisons between Hertz contact loads and contact loads calculated by three- 
dimensional finite element analysis were given in the previous section. In this section, we 
present the effects of the difference in contact loads calculated by two different methods 
on the SIFs. Figures 3.7 through 3.12 compare mode I, II and III SIFs extracted from the 
results of the analyses with Hertz contact conditions and three-dimensional finite clement 
contact analyses. 

SIF values were calculated for the models mentioned in the previous section, 
namely, a pinion with an initial crack, and with cracks corresponding to crack #1 and 


NASA/CR 2003-212529 


29 



crack #2. Two loading conditions were investigated, one corresponding to double tooth 
contact and the other corresponding to single tooth contact, assumed to represent load 
step 3 and load step 1 1 , respectively. 

The difference between SIFs calculated by three-dimensional finite element 
contact analyses and analyses using Hertz contact loads increased as the crack length 
increased due to the enhanced flexibility of the tooth. The change in tooth flexibility as 
the crack grew, resulted in lower contact loads and therefore lower SIFs. 

For the gear/pinion model with an initial crack under double-tooth contact, the 
difference was not expected to be that significant since the initial crack was very small 
and did not introduce significant additional flexibility to the tooth. This behavior is 
shown in Figure 3.7. At crack #1, as shown in Figure 3.8, the change in the SIFs became 
more pronounced. An even more amplified effect of the change in tooth flexibility was 
observed for crack #2, as plotted in Figure 3.9. 
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Figure 3.7: Comparison of mode I, 11 and 111 SIFs calculated by three-dimensional finite 
element contact analysis and finite element analysis with Hertzian contact analysis for a 
pinion with initial crack and at a double tooth contact configuration. 
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Figure 3.8: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 
element contact analysis and finite element analysis with Hertzian contact analysis for a 
pinion with crack #1 and at a double tooth contact configuration. 
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♦ Kllcontact 
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• KHIcontact 
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Figure 3.9: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 
element contact analysis and finite element analysis with Hertz contact loads for a pinion 
with crack #2 and at a double tooth contact configuration. 
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The analyses were repeated for an almost single tooth contact case on all three of 
the crack configurations. The same trends were observed as in the double tooth contact as 
seen in Figures 3.10 through 3.12. However, the difference between the contact analyses 
and analyses with Hertz contact loads were higher for this case. This might be due to the 
fact that the total load applied to the tooth was less than the total load applied in analysis 
with Hertzian contact. As described in section 3.3.1, and shown in Table 3.1, there was at 
least a very little portion carried by the adjacent tooth decreasing the amount of total load 
carried by the cracked tooth. 



■ KIcontact 
□ KIHertz 

♦ Kllcontact 
o KIIHertz 

• KHIcontact 
o Kill Hertz 


Figure 3.10: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 
element contact analysis and finite element analysis with Hertzian contact loads for a 
pinion with initial crack and at a single tooth contact configuration. 

Figures 3.13 through 3.15 show the variation of the ratio of Ku to Kj along the 
crack front for both double-tooth and single-tooth contact. These graphs are plotted using 
the results obtained by FE contact analysis for all three models with initial crack, crack 
#1 and crack #2. 
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Figure 3.1 1: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 
element contact analysis and finite element analysis with Hertzian contact loads for a 
pinion with crack #1 and at a single tooth contact configuration. 



normalized crack length 


Figure 3.12: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 
element contact analysis and finite element analysis with Hertzian contact loads for a 
pinion with crack #2 and at a single tooth contact configuration. 
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Figure 3.13: Ratio of Ku to K/ along the crack front for double-tooth and single-tooth 
contact at initial crack. 



Figure 3.14: Ratio of Ku to K/ along the crack front for double-tooth and single-tooth 
contact at crack #1 . 
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Figure 3.15: Ratio of Ku to Kj along the crack front for double-tooth and single-tooth 
contact at crack #2. 

By performing these contact analyses, we were able to show that the change in the 
stiffness of the tooth due to a growing crack has significant effects on the SIFs and 
therefore on the crack trajectory and fatigue life. If LEFM analyses of section 2 are 
repeated under these new contact loading conditions, the change in the crack trajectory 
and the fatigue life of the pinion, as a result of taking into account the change in the tooth 
flexibility, can be determined. 


4. Summary 

This report summarized new results for predicting crack trajectory and fatigue life 
for a spiral bevel pinion using the Finite Element Method (FEM). The predictions 
presented in this report are based on LEFM theories combined with the FEM, 
incorporating plasticity induced fatigue crack closure and moving loads. 

Previously in this project, three-dimensional computational simulation of crack 
growth in spiral bevel gears was done using BEM [1]. As a result of these analyses, a 
method for predicting three-dimensional, non-proportional, fatigue crack growth 
incorporating moving loads was developed. The BEM studies showed the need for 
improving the efficiency and the accuracy of the computations by employing the state-of- 
the-art computational capabilities. 

In the current work, three-dimensional finite element modeling of the spiral bevel 
pinion was done using OSM/FRANC3D. A three-dimensional FEM mesh of the pinion 
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model was created by a 3D mesher. Moving contact loads on the tooth surface were 
calculated by interpolation using the shape functions on the surfaces of the loaded 
elements. The analyses were carried out using a parallel FEM solver, which calculates 
SIFs using equivalent domain J-integral method. Fatigue life predictions were made 
based on a modified Paris model incorporating crack closure. 

In this report, we show that we can simulate fatigue crack growth in a spiral bevel 
gear more accurately and efficiently using the FEM along with a better representation of 
moving loads. Another very significant improvement was the decrease in solving time of 
the problem by employing parallel PC-clusters. This reduces the computation time to a 
few minutes while substantially increasing the resolution of the solution. 

To obtain a more detailed understanding of the contact between a cracked pinion 
tooth in mesh with an uncracked gear tooth, three-dimensional contact analyses were 
performed on a spiral bevel gear set incorporating a crack. The goal in carrying out these 
analyses was to capture the redistribution of contact loads due to crack growth. Results of 
these analyses showed the expected trend of decreasing tooth loads carried by the cracked 
tooth with increasing crack length. We also showed that this decrease in contact loads 
had an impact on the SIF values and therefore would also affect the crack trajectory and 
fatigue life predictions. 

Predicting crack trajectories more efficiently and accurately is very important for 
gear designers. Development of practical and accurate numerical tools for evaluating the 
gear performance and life enables for new and better gear designs. The work presented 
here provides such a tool capable of predicting crack trajectory and fatigue life in a spiral 
bevel gear. 
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Appendix: Abbreviated sample ABAQUS input file for contact analysis 


♦HEADING 
FRANC 3D Model 

** FRANC3D Version 2.0 Generated Input File 
** NUMBER OF NODES: 13029 
** NUMBER OF ELEMENTS: 58712 
* * 

** NODE DEFINITION 
♦NODE 


1, 3 . 815870/ 0.262049, 1.500000 

2, 3.585600, 0.038081, 1.500000 

3, 3.824640, 0.040620, 1.500000 


13028, 

5.109482, 

0.491277, 

0.734872 

13029, 

5.110899, 

0.478949, 

0.730531 


** END NODE DEFINITION 
* * 

** ELEMENT DEFINITION 
♦ELEMENT, TYPE=C3D4, ELSET=PID1 
1,1904,1855,2100,1903 
♦ELEMENT, TYPE=C3D4, ELSET=PID1 
2,1904,2100,1953,1952 
♦ELEMENT, TYPE=C3D4, ELSET=PID1 
3,125,107,106,846 


♦ELEMENT, TYPE=C3D4, ELSET=PID2 

58711.7717. 12943.7703. 13029 
♦ELEMENT, TYPE=C3D4, ELSET=PID2 

58712. 6078.7717.7703. 13029 
** END ELEMNT DEFINITION 

* * 

♦SOLID SECTION, ELSET=PID1, MATERIAL=MAT1 
♦MATERIAL, NAME=MAT1 
♦ELASTIC, TYPE=ISO 
3e+07 , 0.3 
♦DENSITY 
7800 

♦SOLID SECTION, ELSET=PID2, MATERIAL=MAT2 
♦MATERIAL, NAME=MAT2 
♦ELASTIC, TYPE=ISO 
3e+07 , 0.3 
♦DENSITY 
7800 
* * 

♦♦Definition of contact for tooth 1 
* * 

♦ELSET, ELSET=MASTER_1 
26454, 

26457, 

26458, 
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57531, 

58366, 

58367, 

*ELSET, ELSET=SLAVE_1 
1 / 

2 , 

5 , 


22208, 

24437, 

24438, 

♦SURFACE DEFINITION, NAME=MASTER_1 

26454. 54 

26457. 54 

26458. 54 


57531. 51 

58366. 51 

58367. 51 

♦SURFACE DEFINITION, NAME=SLAVE_1 

1, S2 

2 , S4 
5, S4 


22208. 51 

24437. 51 

24438. 51 

♦SURFACE INTERACTION, NAME=NOFRIC 
♦FRICTION 
0 

♦CONTACT PAIR, INTERACTION=NOFRIC 
MASTER_1 , SLAVE_1 
* * 

** Definition of contact for tooth 2 
* * 

♦ELSET, ELSET=MASTER_2 
26455, 

26456, 

26459, 


57320, 

57753, 

57754, 

♦ELSET, ELSET=SLAVE_2 

3 / 

4, 

9 , 
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24532, 

24540, 

24548, 

♦SURFACE DEFINITION, NAME=MAS TER_2 

26455. 54 

26456. 54 

26459. 54 


57320. 51 

57753. 51 

57754. 51 

♦SURFACE DEFINITION, NAME=SLAVE_2 

3, S2 

4 , S4 
9, S4 


24532. 51 

24540. 51 

24548. 51 

♦SURFACE INTERACTION, NAME=NOFRICl 
♦FRICTION 
0 

♦CONTACT PAIR, INTERACTION=NOFRICl 
MAS TER_2 , S LAVE_2 
* * 

** Application of displacement to ensure contact 
** between gear and pinion 
* * 

♦STEP, NLGEOM, INC=100 

♦♦FRANC3D load case 1 

♦STATIC 

1 . 0 , 1.0 

♦BOUNDARY 

2209, 3, 3, 0 

2211.3.3, 0 

2213.3.3, 0 

2215.3.3, 0 

2217.3.3, 0 

2219.3.3, 0 

2221.3.3, 0 
2223, 3, 3, 0 

2232.2.2, -0.01 

2247.3.3, 0 


33.1.1.0 

33.2.2.0 

33.3.3.0 
2322, 1, 1,0 
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2322.2.2.0 

2322.3.3.0 

2415.1.1.0 

2415.2.2.0 

2415.3.3.0 
*END STEP 
* * 

** Application of full design torque 
* * 

*STEP, NLGEOM, INC=100 
**FRANC3D load case 1 
♦STATIC 
1 . 0 , 1.0 

♦BOUNDARY, OP=NEW 
2209, 3, 3, 0 

2211.3.3, 0 

2213.3.3, 0 

2215.3.3, 0 

2217.3.3, 0 

2219.3.3, 0 

2221.3, 3, 0 

2223.3.3, 0 

2247.3.3, 0 


33.1.1.0 

33.2.2.0 

33.3.3.0 
2322, 1, 1,0 

2322.2.2.0 
2322, 3, 3,0 

2415.1.1.0 

2415.2.2.0 

2415.3.3.0 
♦CLOAD 

2799.2, -1000 
♦END STEP 

♦STEP, NLGEOM, INC=100 

♦♦FRANC3D load case 1 

♦STATIC 

1 . 0 , 1.0 

♦CLOAD 

2799.2, -3230 
♦FILE FORMAT, ASCII 
♦NODE FILE, FREQUENCY=100 
U 

♦EL FILE, FREQUENCY=100, POSITION=NODES 
S, E, PE 

♦RESTART, WRITE, FREQUENCY=100 
♦NODE PRINT, FREQUENCY=100 
U 

♦EL PRINT, FREQUENCY=100, POSITION=NODES 

♦CONTACT PRINT, FREQUENCY=100 

CFN, CARE A, XN 

♦END STEP 

♦♦END OF INPUT 
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